/*-------------------------------------------------------------------------*/
/**
  @file		randnum.c
  @author	D E Knuth. 
  @date	july 2004
  @version	$Revision: 2.2 $
  @brief	portable random number generator
  
  http://www-cs-faculty.stanford.edu/~knuth/taocp.html

*/
/*--------------------------------------------------------------------------*/



/*controlling macro*/
#include "plib.h"

#define KK 100                     /* the int lag */
#define LL  37                     /* the short lag */
#define MM (1L<<30)                 /* the modulus */
#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */

int ran_x[KK];                    /* the generator state */

#ifdef __STDC__
void ran_array(int aa[],int n)
#else
void ran_array(aa,n)    /* put n new random numbers in aa */
  int *aa;   /* destination */
  int n;      /* array length (must be at least KK) */
#endif
{
  register int i,j;
  for (j=0;j<KK;j++) aa[j]=ran_x[j];
  for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
  for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
  for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
}

/* the following routines are from exercise 3.6--15 */
/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */

#define QUALITY 1009 /* recommended quality level for high-res use */
int ran_arr_buf[QUALITY];
int ran_arr_sentinel=-1;
int *ran_arr_ptr=&ran_arr_sentinel; /* the next random number, or -1 */

#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
int ran_arr_cycle()
{
  ran_array(ran_arr_buf,QUALITY);
  ran_arr_buf[100]=-1;
  ran_arr_ptr=ran_arr_buf+1;
  return ran_arr_buf[0];
}

#define TT  70   /* guaranteed separation between streams */
#define is_odd(x)  ((x)&1)          /* units bit of x */

#ifdef __STDC__
void ran_start(int seed)
#else
void ran_start(seed)    /* do this before using ran_array */
  int seed;            /* selector for different streams */
#endif
{
  register int t,j;
  int x[KK+KK-1];              /* the preparation buffer */
  register int ss=(seed+2)&(MM-2);
  for (j=0;j<KK;j++) {
    x[j]=ss;                      /* bootstrap the buffer */
    ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
  }
  x[1]++;              /* make x[1] (and only x[1]) odd */
  for (ss=seed&(MM-1),t=TT-1; t; ) {       
    for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
    for (j=KK+KK-2;j>=KK;j--)
      x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
      x[j-KK]=mod_diff(x[j-KK],x[j]);
    if (is_odd(ss)) {              /* "multiply by z" */
      for (j=KK;j>0;j--)  x[j]=x[j-1];
      x[0]=x[KK];            /* shift the buffer cyclically */
      x[LL]=mod_diff(x[LL],x[KK]);
    }
    if (ss) ss>>=1; else t--;
  }
  for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
  for (;j<KK;j++) ran_x[j-LL]=x[j];
  for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
  ran_arr_ptr=&ran_arr_sentinel;
}

/*swap: swap the values of two pointers.
  Pre: p1 and p2 are valid pointers.
  Post: the values of p1 and p2 have already been swapped.
*/
void swap(int *p1, int *p2)
{
	int temp;
	temp = *p1;
	*p1 = *p2;
	*p2 = temp;

}


/*randnum: Generate consecutive random numbers in a given distribution.
  Pre: the distribution is given by [0, range], where range >=0, a[] is the array to store the result.
  Post: consecutive random numbers have been already generated in the array a[].
  Use: swap().
*/
void randnum(int range, int *a)
{/*cgseed = -1 means there is no need to have random seed.*/
	int i;
	int tmp;
	
	for (i=0; i <range; i++)
		a[i] = i;

	for (i=0; i < range; i++){
		tmp = (int)fmod(abs(ran_arr_next()), range);
		swap (&a[i], &a[tmp]);   
	}
}
